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Abstract 



The vibration suppression effectiveness of a flexible in-line marine cable vibration 
absorber is studied. The transfer matrix method is used to build various numerical models 
of vibration absorbers in marine cable systems. The models determine cable system 
natural frequencies, mode shapes and modal damping ratios. 

The introduction of absorber damping is shown to result in complex roots to the 
modal characteristic equations. A computer complex root solver is used to solve for the 
complex roots of the characteristic equations, resulting in complex system natural 
frequencies. The significance of complex natural frequencies is explained. Complex 
natural frequencies are used to calculate modal damping ratios. 

The models demonstrate that absorber effectiveness is heavily dependent on 
absorber location, absorber mass and absorber length. Parametric variation is used to 
achieve maximum effectiveness of the flexible in-line absorber. Even under optimum 
conditions, it is shown that the absorber provides insufficient damping to reduce vortex- 
induced vibrations in water. 

The same transfer matrix method is used to evaluate the effectiveness of a mass- 
spring-dashpot type absorber in a marine cable system. This type of absorber is shown to 
produce adequate damping to reduce vortex-induced vibrations in water. The transfer 
matrix method used in this thesis is validated by analyzing the same system using an 
approach by Den Hartog [1]. The transfer matrix approach combined with complex root 
solving capability is shown to provide an effective analysis method for marine cable 
systems. 
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Chapter 1 
Introduction 



Long slender marine structures such as cables and pipelines are subject to a 
condition called vortex induced vibration. This occurs when the cylindrical structure is 
subjected to a cross flow current as shown in figure 1-1: 




Figure 1-1: Vortex Shedding 

At a critical value of flow speed, U, a periodic shedding of vortices occurs on the upper 
and lower surfaces of the cylinder. This behavior is known as a Von Karman street. As 
each vortex is shed into the wake, a localized pressure drop occurs. The result is a 
harmonically varying vertical (lift) force on the cylinder. The frequency at which vortex 
shedding occurs is related to cylinder diameter and flow velocity by a dimensionless 
parameter called the Strouhal number; 

f D 

St= y (1-1) 



where: f = vortex shedding frequency (Hz) 

D = cylinder diameter (m) 

U = flow speed (m/sec) 
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The Strouhal number is approximately 0.2 for a wide range of Reynolds numbers. 
Vortex-induced vibrations are especially severe when the vortex shedding frequency is 
near a resonance frequency of the structure. Cross-flow displacement of the cylinder at 
this “lock-in” frequency may be on the order of one cylinder diameter. Such vibrations 
can fatigue marine structures and cause eventual failure. 

This problem motivates the investigation into a dynamic absorber for long 
cylindrical marine structures. Before analyzing the effect of various damping devices, it is 
worth while to determine the magnitude of damping which must be introduced into a 
marine structure to cause significant reduction in vortex induced vibrations. In figure 1-2 
below, Griffin [13] shows how cross flow displacement varies with a dimensionless 
parameter called the reduced damping. So. Reduced damping is the ratio of damping 
force to fluid exciting force. 




Figure 1-2: Vibration Displacement vs. Reduced Dani|)ing 
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where; 



Sg = 8 71" St^ ( Ills / Pf D^) s reduced damping 



(1-2) 



St = Strouhal number 

ms = linear density of cable (kg/m) 

pf = fluid density (kg/m'’) 

D = cable diameter (m) 

The reason displacement amplitude approaches a plateau value at low reduced damping 
(values below Sg about 0.2) is due to the self-limiting nature of vortex induced vibrations. 
At values of displacement amplitude greater than approximately one cylinder diameter, 
the extreme motion of the cylinder begins to disturb the formation of the Von Karman 
street. Stuctural damping is not the limiting factor at low values of reduced damping. The 
reduced damping is valuable in determining the amount of damping which must be added 
to the structure to have an effect on displacement amplitude. It would not be worthwhile 
to add damping such that reduced damping would not exceed the “knee in the curve” 
value of 0.2 illustrated in figure 1-2. This is also a very convenient approach since the 
damping ratio used to calculate reduced damping is the damping which one would 
measure in an in vacuo transient decay test. Thus, various absorber configurations can be 
evaluated without the need to consider hydrodynamic effects 

This thesis will study various vibration absorber configurations for marine cables. 
Most of the analysis will focus on the type of absorber shown in figure 1-3 on the 
following page; 
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cable 



damper 

1/ 




Rigid masses 



Figure 1-3: Cable Vibration Absorber 

This absorber was first proposed by Vandiver and Li [10] . Energy dissipation occurs by 
way of relative angular motion at the pivot between the two rigid masses. The damping 
device is integral to the pivot but it is shown schematically in figure 1-3. Butler [7] 
showed this type of damping device to be possible for marine cable applications. It was 
seen as a very versatile configuration for marine cables since such an absorber could be 
clamped on at a desired location while deploying a cable. Or, if the absorber were 
permanently installed, its size and orientation might permit the absorber to be reeled onto 
a cable drum when the cable is not deployed. 

Chapter 2 describes the transfer matrix approach for analyzing multi-component 
mechanical systems M’ith finite boundaries. It is a very convenient method for studying 
the dynamic behavior of marine cable systems with attached damping devices. The 
analysis utilizes the absorber transfer matrix derived by Li [5]. 

In chapter 3, the transfer matrix method is used to determine natural frequencies, 
mode shapes and damping ratios of the cable / absorber system. A method is formulated 
for determining the optimum value of damping constant. 

In chapter 4 , system parameters such as absorber length and mass are varied to 
study their effect on damping performance. Also, use of multiple absorbers and the effect 
of absorber location relative to vibration nodes and antinodes was studied. Use of a more 
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conventional absorber type is studied by merging the transfer matrix approach with a more 



classical approach (Den Hartog [1]) with some interesting results. 

Chapter 5 concludes the study and suggests areas for further research. 
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Chapter 2 



The Transfer Matrix Method 

The transfer matrix method allows for a convenient matrix representation of the 
marine cable system. In 1963, Pestel and Leckie [2] recognized that the introduction of 
numerical computing made matrix methods a powerful tool for studying mechanical 
system dynamics. In marine cable systems, features such as rigid bodies and discrete 
damping devices are well suited to matrix modeling techniques. Today’s computing speed 
allows complicated marine cable systems to be modeled with relative ease. 

2.1 The Marine Cable System 

2.1.1 Assumptions 

It is important here to state the assumptions employed in representing the cable 
system in this way. These assumptions are those of the classic linear string wave equation. 
They were used by Li [5] in developing the transfer matrices used in this analysis; 

1) the cable systems are taut and undergo only small deflections 

2) cable bending stiffness is neglected 

3) the masses in the cable system are rigid 

4) the weight of these masses is small compared to cable tension 

5) the dynamic variation in tension in the cable is small and neglected 

6) only planar motion is considered (2 dimensional) 
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2.1.2 Transfer Matrix Representation of a Marine Cable System 

Transfer matrices are used to relate the physical state of one location in a 
mechanical system to that at another location in the system. The state at any location in 
the system is represented by the state vector. In a cable system, the state at any location 
can be represented by the state vector, z: 

bl= 

where y is the Fourier transform of cable transverse displacement and y' is the Fourier 
transform of cable slope. The product of cable tension T and'y’ is a value equal to the 
transverse force caused by tension at the location of the state vector. Thus, the cable 
system state vector contains a displacement value and a force value. If assumption 2 in 
section 2.1.1 were not valid and cable bending stiffness could not be neglected, the state 
vector would require four elements: displacement, slope, shear force and bending moment. 
The use of Fourier transforms is useful since the types of motion studied here will limited 
to free and forced harmonic vibrations. 

Within a mechanical system, the state vector at any location can be obtained by 
multiplying the state vector at some other location in the system by the transfer matrix of 
the component(s) between the two locations. In a marine cable system, the transfer 
matrices are 2 by 2 matrices represented as follows: 



y(x) 

TT(x) 



(2-1) 



The elements of the transfer matrices are derived using the equations of motion. Since 
they relate the Fourier transforms of displacement and force, the transfer matrix elements 
are functions of circular frequency, o. The transfer matrix of a component is general in 
that it does not change when the component is installed in other systems. Cable system 
components modeled here using transfer matrices include cable segments, lumped and 
distributed masses and vibration absorbers. 

Figure 2-1 illustrates how transfer matrices are used to relate the state vectors at 
individual locations in the cable system; 



/ 



/ 



ty 



[T2] 



[T.] 



Xl 



X2 



absorber 



[T:0 



\ 



X', 



[zi] 



O 



[Z2] 

© 



[Z3] 

© 



X 4 \ 

\ 

[Z4] 



© 



[Z2] = [Ti] [z,] 


(2-3) 


[Zs] = [Ti] [T 2 ] [Zi] 


(2-4) 


[Zs] = [T 2 ] [Z2] 


(2-5) 


[z(x)] = [Ti(x)] [z,] ; X <X2 


(2-6) 



Figure 2-1: Transfer Matri.x Representation of Cable / AI)sorl)er Sssteni 

The transfer matrix for the uniform cable section, [Ti], comes directly from the 

solution to the linear string wave equation and is given by; 

cos(k:/) l/kTsin(k/) 

-kTsin(k/) cos(k/) (- ^) 
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where: 



k = wave number = co / c (m'‘) 
c = wave propagation speed = (T / p)' * (m/s) 

T = cable tension (N) 
p = cable linear mass density (kg/m) 

/ = cable length from Xi to X 2 (m) 

For locations on the cable to the left of X 2 , substitution of x for / in equation 2-7 will yield 
the transfer matrix [Ti(x)] as in equation 2-6. 

The transfer matrix for the absorber, [T 2 ], is given in Appendix B and in the thesis 
by Li [5], Appendix B also contains the transfer matrices for various vibration absorbers 
and discrete masses which may be located within a cable system. 

2.1.3 The Overall Transfer Matrix 

The product of the three transfer matrices, [Ti] [T 2 ] [T?], equals the overall 
transfer matrix of the system. For the system shown in figure 2-1, the boundaries are 
walls. For a wall, cable displacement would be zero and cable force would be unknown. 
Thus, the state vectors at the boundaries can be represented by: 




The overall system transfer matrix can be used to relate the two boundary state vectors: 
[T0V..11] = [Ti] [T2] [T] (2-9) 

[Z4] = [Toverall] [Zl] (2-10) 
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(2-11) 



'0 




“til ti:' 


'0 


.Ty'(x 4 )_ 




t 2 i t:: 


TT(xi)_ 



The zero displacement at the boundaries results in; 

(tiix 0 )+(ti:xTy’(xi)) = 0 ( 2 - 12 ) 

=^>ti:=0 (2-13) 

This is the characteristic equation for the system. Since transfer matrix elements are 
functions of frequency, the roots of t^ of [Toveraii] are the natural frequencies of the 
system. 

The mode shapes of the system can then be plotted by setting a unity force at the wall: 

'o' 



[z,] = 



(2-14) 



It follows that: 

[z(x)] = [Ti(x)] [zi] ; x<X 2 (2-15) 

-[Ti][T2][T3(x)][zi];x>X3 (2-16) 

Note that [T. 3 (x)] is expressed by substituting the quantity (x - X 5 ) for / into equation 2-7. 
This coordinate transformation is required because transfer matrix expressions normally 
reference a coordinate system which is located on the component that the transfer matrix 
represents. 

In general, 



[z(x)] = [T(x)] [z,] = 



tii(x) ti:(x) 


'o' 


t:i(x) t::(x) 


_ 1 _ 



(2-17) 



where [T(x)J is the transfer matrix product of the system components between the left wall 
and X. 
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The first element of [z(x)] is therefore; 

y(x) = ti;(x) (2-18) 

Equation 2-18 specifies the displacement mode shape of the system. Note that ti 2 (x) in 
equation 2-18 is an element of the transfer matrix between the left wall and an arbitrary 
point, X, and is not the same as tn of ToverMi in equation 2-13. 

In chapter 3, this approach will be applied to a specific marine cable / absorber 

system. 
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Chapter 3 



Vibration Absorber Analysis 

As mentioned in the introduction, the vibration absorber which will be the primary 
focus of this analysis is that shown below in figure 3-1 : 




Figure 3-1: In Line Vibration Absorber 

The transfer matrix expression for this absorber was derived by Li [5] and is quite lengthy, 
requiring a page and a half to state (see Appendix B). This is because the transfer matrix 
must account for both the translational and rotational motion of the masses, damping and 
stiffness. MATLAB was used to compute the values of this transfer matrix and was used 
for all other modeling and calculation in this thesis. Appendix A contains the MATLAB 
programs STRING 1 and STRING2 which calculate the transfer matrices of the absorber 
and cable segments, respectively. 

3.1 Modeling the Vibration Absorber in a Cable System 

Since the nature of this analysis was numerical, a baseline absorber system of the 
type shown in figure 3-1 was modeled. System parameters such as cable length, cable 
diameter and absorber mass were set to values which could realistically be expected in a 
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marine system. Optimum absorber performance can be obtained by systematically varying 
these parameters. Baseline system parameters are as follows: 



Total system length: 


50.8 m 


Cable tension, T: 


5000 N 


Cable density, p: 


1 kg/m 


Absorber link mass, M: 


9 kg 


Absorber link moment of inertia, I: 


0.12 kg m^ 


a and b dimensions: 


0.03 m 


c, d, e, f dimensions: 


0.2 m 


Total absorber length, Labs' 


0.8 m 


Stiffness, K: 


lON/m 


Damping, R: 


varies N/(ni/s) 



Note that the stiffness, K, is relatively insignificant. Cable tension is the dominant factor in 
providing a restoring force for the absorber 

Absorber analysis consisted of the following steps: 

1) Obtain natural frequencies 

2) Determine damping ratio 

3) Determine displacement mode shapes 

4) Vary system parameters to optimize damping ratio 
This process is detailed in the following sections. 
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3.2 The Cable / Absorber System With no Damping 



It is valuable to first analyze the system and demonstrate the analysis process with 
the damping value of the absorber set to zero. This keeps all calculations in the real 
domain and provides for a better understanding of the introduction of damping. 

3.2.1 System Natural Frequencies 

Chapter 2 detailed the process for determining the characteristic equation and 
natural frequencies using transfer matrices. The overall transfer matrix of the system, 
[Toveraii], is computed as follows; 

[T„v.„ii] = [T,] [To] [Tj] (3-1) 

where; 

[Ti] and [T?] are the transfer matrices of the left and right cable segments 
as in equation 2-7 

[T 2 ] is the transfer matrix of the absorber as in Appendix B 
Due to the fixed wall boundary conditions, the system characteristic equation is 
ti2=0 (3-2) 

where ti 2 is the first row, second column element of Toveraii- 

A MATLAB script file, STRING4, was used to generate figure 3-2 on the 
following page, a plot of ti 2 versus frequency, to. The zeroes of ti 2 are the natural 
frequencies of the system. 
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T(1,2) CHARACTERISTIC FOR SYSTEM WITH ZERO DAMPING 




Figure 3-2: ti 2 Characteristic of Uiulani|)ed System 

3.2.2 System Mode Shapes 

Once the system natural frequencies are known, transfer matrix multiplication can 
be used to determine the displacement and force mode shapes. This is accomplished by a 
MATLAB program called STRING6 which takes the system frequency as an input. The 
state vector at the left wall is set as follows 



0 



Zieft 



I 



( 3 - 3 ) 



where the first element indicates the zero displacement at the wall and the second element 
is set at one. This results in a mode shape which is normalized to a unity transverse force 
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at the system boundaries. The displacement mode shape is calculated along the length of 
the system as described in Chapter 2: 
from left wall to absorber; 



[z(x)] = [Ti(x)jzieft] 



til tl2 


'o' 




’y(x) 


t21 t22 


_ 1 _ 




_-Iy'(x)_ 



y(x) = ti:(x). 



where ti 2 (x) is an element of [Ti(x)] 
from absorber to right wall: 
y(x) = ti:(x) , 



(3-4) 



(3-5) 



(3-6) 



where ti 2 (x) is an element of the matrix product [Ti] [T 2 ] [T:,(x)]. 
Similarly, the transverse force mode shape would equal the t 22 element of the transfer 
matrix as a function of distance from the left wall. The resulting displacement mode 
shapes for the first three vibration modes are shown in figure 3-3 on the following page. 
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Figure 3-3: Undamped System Displacement Mode Shapes 
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3.3 Introduction of Absorber Damping 



The introduction of absorber damping takes the transfer matrix calculations into 
the complex domain. In addition to natural frequency and mode shape, damping ratio can 
now be determined. Damping ratio is the parameter used in determining absorber 
effectiveness. Damping ratio is used to compute the reduced damping in figure 1-2 and to 
determine if sufficient damping exists to reduce vortex-induced vibrations. 

3.3.1 System Natural Frequencies 

The same approach which was used to compute undamped natural frequencies in 
section 3.2.1 is used here to compute damped natural frequencies. As shown in section 
3.2.1, the characteristic equation of the system is: 

ti2=0 (3-7) 

where tn is an element of the overall system transfer matrix. Figure 3-4 on the following 
page shows the magnitude of ti 2 as a function of frequency, co. Note that when damping 
is present, ti 2 is a complex quantity. The figure also shows ti 2 for the undamped system 
in the same frequency range 

The damped curve is quite different in appearance from the undamped case. For 
alternating modes, the damped value of tn fails to reach a zero for real values of natural 
frequencies. This behavior is significant and will later be shown to have a very physical 
explanation. 
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T(1.2) T(1.2) 




R=8000 N/(m/s) 




w(rad/sec) 



Figure 3-4: Magnitii(Ie(ti 2 ) Versus Frec|iieiicy for Damped ss. Undamped System 
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In order to obtain all the zeroes of ti 2 , complex values of frequency must be considered. 
This requires the use of a computer complex root solver, MULLER (Appendix A), 
developed by Conte and de Boor [14] and converted into MATLAB format by Wilson and 
Roberts. Before proceeding further with the results of the complex root search, the 
meaning of a complex natural frequency must first be discussed. 

The harmonic motion at a point in the cable system can be described by; 
y(x) = Ae’“', (3-8) 

where A is displacement amplitude. 

For imaginary natural frequencies; 

CO COreal ^COimag (j-9) 

Substituting; 

y(x) = (3-10) 

— ^ gi(wreal)t g-(o)iiuag)t (3-11) 

= A[COS(0realt) + isin(cOrealt)] (3-12) 

The real part of which is; 

= A cos(o)„.,0 e-'"'”-'*' (3-13) 

So, we see that for free vibrations the damped system vibrates at co,eai with a decay 
envelope The ratio of the real and imaginary parts of the natural frequency is the 

modal damping ratio; 

C — COimag/ COreal (-^*1^) 

Now that the significance of a complex natural frequency is understood, we can 
return to the analysis of the cable / absorber system. The complex root solver, MULLER, 
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is iterative and requires a starting guess for the complex frequency root. The undamped 
natural frequencies were used as initial guesses when running MULLER. With the 
assistance of Dr. Rama Rao, MIT, a routine called MULRUN (Appendix A) was 
developed which successively runs MULLER for many values of damping constant to 
obtain an output which demonstrates the optimum value of damping. Figure 3-5 below 
shows the output of MULRUN for the first vibration mode of the baseline system. 




Figure 3-5: Demoiisti atioii of Optimum Damping Value 
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0 



The MULRUN output shows that as damping constant, R, is varied, there is an 
intermediate damping value which optimizes (maximizes) damping ratio. It was theorized 
that the reason damping ratio decreases at high values of damping constant is that relative 
motion between the absorber masses decreases above the optimum value. This was 
confirmed by creating a transfer matrix model of a uniform bar (see Appendix B) in place 

of the vibration absorber in the baseline system as shown in figure 3-6: 

damper 





cable 



Uniform bar 




Figure 3-6; Representing High Stiffness Absorber as Uniform Bar 
The uniform bar had the same dimensions and mass as those of the absorber. It was 

shown that the cable system with the uniform bar at its center had the same natural 
frequency and mode shape as the absorber with very high damping constant. So at very 
high damping, relative rotational motion between the absorber links is restricted and the 
absorber is unable to dissipate energy. The absorber acts as a uniform bar. 

3.3.2 System Mode Shapes With Damping 

Once the system natural frequencies are known, transfer matrix multiplication can 
be used to determine the displacement and force mode shapes using the same approach as 
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for the undamped case. This is accomplished by a MATLAB program called STRING6 
which takes system frequency as an input. For the damped case, complex natural 
frequencies are input into STRING6. The state vector at the left wall is again: 




The displacement mode shape is calculated along the length of the system as in section 
3.2.2. The resulting mode shapes for the first three vibration modes are shown in figure 
3-7. 






Figure 3-7: Damped System Displacement Mode Shapes 
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Inspection of these mode shapes reveals that for the even mode, the absorber is located at 
a modal antinode. There is no relative motion between the absorber links. The absorber 
simply rotates about its center as if it were a uniform bar. This is confirmed further by 
using MULRUN for the same modes to determine optimum damping ratio. In figure 3-8 
below, it is seen that for the even vibration mode, damping ratio is insignificant even 
though damping is present. Relative motion between the absorber links is required if the 
absorber is to dissipate energy. 



First Mode 






Figure 3-8: Damping Effect in First Three V'ihratioii Modes 
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So, from initial analysis the following observations are made: 

1) An optimum value of damping constant is one which adds sufficient 
damping while not overly restricting motion between the two absorber 
links. 

2) Absorber effectiveness is very dependent on absorber location. 
Absorbers located at modal antinodes provide the most energy dissipation. 

These observations are consistent with the behavior of more conventional types of 
vibration absorbers, such as those installed in discrete (vice continuous) mechanical 
systems. 

In chapter 4, parameters such as absorber mass, absorber length and number of 
absorbers will be varied to study their effect on absorber performance. 
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Chapter 4 



Parameters Affecting Absorber Performance 

In previous chapters, it has been demonstrated how natural frequencies, mode 
shapes and damping ratios can be obtained for a vibration absorber in a marine cable 
system. It was shown how the damping constant, R, can be varied to arrive at an 
optimum damping ratio. The tools are now in place to study the effects of other 
parametric variations on absorber performance. Improvements in absorber performance 
will be sought by varying the following: 

1) Absorber length 

2) Absorber mass 

3) Number of absorbers installed in a single cable system 

4) Absorber location 

4.1 Required Damping 

Before proceeding with further attempts to optimize damping, it would be helpful 
to know what value of damping is required in the system to reduce vorte.x induced 
vibrations. Recall from figure 1-2 that there is a value of reduced damping below which 
structural damping provides no assistance in reducing vortex induced vibrations. 

Structures with reduced damping in this plateau range will have displacement amplitudes 
which are restricted by limit-cycle reductions to the lift forces. Thus, there is no gain in 
adding small amounts of damping which will not raise reduced damping above that value 
required to reduce the response to less than that resulting from limit cycle behavior. From 
figure 1-2, this threshold value of reduced damping is as follows: 



Sg>0,2 



where; 



Sg = 8;c- St'C(ms/pfD‘) 



(4-1) 



For the baseline system in water (see section 3.1), this results in a threshold (required) 
value of damping ratio; 



The highest value of damping ratio achieved thus far (see figure 3-8) is approximately 
Q = .0012. So, we see that if the absorber being studied is to be of use in reducing vortex 
induced vibrations in water, further optimization is required. 

4.2 Variation of System Parameters 

4.2.1 Absorber Length 

In figure 3-7, mode shapes for the first three modes were shown. These low 
modes have very long wavelengths, resulting in the situation shown in figure 4-l(a) on the 
following page. When absorber length is smaii compared to wavelength, relative angular 
motion between the absorber links is very small, even in odd numbered vibration modes. 
Since energy dissipation relies on this relative motion, short absorber links result in poor 
absorber performance. As absorber length is increased, relative motion between the 
absorber links is increased, and energy dissipation increases as shown in figure 4- 1(b). A 
key dimensionless parameter of use here would be the ratio of absorber length to the half 
wave length of a given mode. This parameter is designated as the length ratio: 



C>.04 



(4-2) 




(4-3) 
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Figure 4-1; Effect of Absorber Length on Dam|>er Motion 
Variation of absorber length was accomplished in an indirect way. Instead of varying the 

absorber link length, the effect of absorber length was studied by varying X, or simply by 
comparing the damping ratios of different modes of the baseline system. Starting with 
mode 1, optimum damping values were determined for each mode. As discussed in 
section 3.3.2, even numbered vibration modes e.xhibited no damping effects. The best 
absorber performance (as determined by ma.ximum damping ratio) was achieved in the 
29th mode. 

In this mode: 

X = 3.50 m 

Labsorbei -8 m 

LR= .457 

In modes beyond mode 29, higher length ratio resulted in more restricted absorber motion 
and reduced damping ratios. 



Figure 4-2 shows the damping value optimization and mode shape for mode 29. 



DAMPING RATIO VS. DAMPING CONSTANT. R 





Fi^re 4-2: Damiiing Optimization anti Mode Sha|)c for Mode 29 
It can be seen that with an optimum value of damping constant, damping ratio is 
approximately 1.7%. This falls well short of the value of 4% required to reduce vortex 
induced vibration amplitude in water. 

4.2.2 Absorber Mass 

In further efforts to increase the amount of damping which the absorber adds to 
the system, the mass of the absorber links was varied. A useful dimensionless parameter 
for this analysis was the mass ratio 

hTR. m^t^sorber 7 (I-absorber Pcable ) (^“^) 

This is the ratio of the absorber mass to the mass of the cable segment displaced by the 
absorber. 
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The variation of mass ratio was performed for the 29th mode of vibration of the 
cable system. This mode was selected because it this was the mode for which optimum 
absorber length was shown . Thus the result of this variation represents an attempted 
optimization of both absorber mass and length. There is no guarantee that optimum length 
ratio remains at 0.457 when mass ratio is varied. Figure 4-3 shows how optimum 
damping ratio varies with absorber mass: 



Mass Ratio - 16 





Mass Ratio = 22.5 




Fij^ure 4-3: Affect of Muss Ratio on Mode 29 Dam pi a 
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Figure 4-3 shows that increasing absorber mass ratio tends to increase optimum 
damping ratio. Also, lower values of damping constant are required to achieve the 
optimum damping value when a higher mass ratio is used. This figure demonstrates that 
the baseline system (MR=22.5) has the optimum mass ratio because it results in the 
highest system damping ratio. When mass ratio was increased beyond 22.5, the system 
shifted to an even numbered vibration mode where damping is insignificant. 

If this type of analysis were to be performed on a real marine system, one would 
have to exercise care in changing absorber mass. Changing absorber mass changes the 
mass of the vibrating system and thus changes system natural frequencies. The marine 
system would be designed with some range of current velocity and corresponding vortex 
shedding frequencies in mind. If one were not careflil, one could change system mass such 
that the absorber was optimized at some frequency outside the range of concern. 

4.2.3 Behavior at Realistic Values of Absorber Mass 

Admittedly, a mass ratio of MR=22.5 corresponds to a very large absorber mass, 
equaling 36% of the total cable mass. Absorbers of this extreme size were considered 
only in an attempt to achieve the required damping level of section 4.1. A more typical 
absorber mass might be no more than 10% of the cable mass. This would correspond to a 
mass ratio of MR=6.25. In this section, an absorber of this size was analyzed over many 
vibration modes and corresponding length ratios. The results are shown in figure 4-4 and 
4-5 on the following page. 
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OPTIMIZED DAMPING RATIO VS. MODE 
NUMBER AND LENGTH RATIO 




Length Ratio and Corresponding Mode Number 



Figure 4-4: Optimized Damping Ratio \s. Mode Number and Length Ratio for Mass Ratio = 6.25 

OPTIMUM DAMPING CONSTANT VS. MODE 
NUMBER AND LENGTH RATIO 




Figure 4-5: Ojitimized Damping Constant \s. Mode Number and 
Length Ratio f(»r Mass Ratio = 6.25 

Figures 4-4 shows that an optimum damping ratio is achieved in mode 53 with a length 
ratio of LR=0,83. Figure 4-5 shows how required damping decreases dramatically with 
increasing length ratio. This behavior was discussed in section 4.2.1. 
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4.3 Number of Absorbers 



Up to this point, variation of system parameters has not resulted in achieving the 
required damping ratio of Section 4.1. It would be very convenient if one could 
multiplicatively increase damping ratio by simply increasing the number of absorbers 
installed. 

This addition of absorbers was performed for mode 3. A transfer matrix model 
was constructed for a cable system with vibration absorbers at its mode 3 antinodes. An 
absorber mass ratio of MR=6.25 was used. All other system parameters were those of the 
section 3.1 baseline system. The MATLAB programs which model this system are located 
in directory MULTABS in Appendix A. These programs are functionally identical to 
those used to evaluate the system with one absorber. 

Figure 4-4 on the following page shows the mode 3 mode shapes for the system 
with one absorber and the system with three absorbers. For both cases, a mass ratio of 
MR=6.25 and a length ratio of LR=.047 was used. The maximum damping value 
achieved in each case is indicated 
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^ ^q- 3 3rd MODE SHAPE. ONE ABSORBER AT CENTER 




3rd MODE SHAPE, DAMPER AT EACH ANTINODE 




Fij^ure 4-6: Comparison of Single Absorber to Multiple Absorbers 

Inspection of figure 4-6 reveals that optimized damping ratio roughly doubled 
when the number of absorbers was tripled. Thus, damping ratio is not multiplicative with 
the number of absorbers installed in the system. One should also observe from figure 4-4 
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that mode 3 natural frequency decreases considerably with the additional mass loading of 
the extra absorbers. 

4.4 Absorber Location for Operation in Multiple Modes 

In previous sections, absorber performance has been looked at one mode at a time. 
Parameters such as damping constant and absorber mass have been varied and the effect 
on modal damping ratio has been observed. This has improved understanding of the 
system but it is not the way in which a marine vibration absorber system would be 
designed. In a system subject to vortex induced vibrations, some range of potential 
currents would exist. This range of currents would be directly correlated to an excitation 
bandwidth by the Strouhal relationship of equation 1-1. So, an absorber design problem 
would be concerned with multiple modes. This section will illustrate how the proper 
location of a single absorber can optimize performance in multiple vibration modes. 

In section 3.3.2, it was observed that the most effective location for an absorber 
was at a modal antinode. For mode 3, these locations would be L/6, L/2, and 5L/6. Now 
suppose we wanted to design an absorber system which would be effective for modes 3 
and 4. The mode 4 antinodes are at L/8, 3L/8, 5L/8. The two modes share no common 
antinodes. In fact, the L/2 antinode for mode 3 is a mode 4 node. So, an absorber located 
at L/4 offers no absorber effect for mode 4. 

A compromise was arrived at by constructing a model of a system with an 
absorber at 7L/16. This location is half way between the 3L/8 mode 4 antinode and the 
L/2 mode 3 antinode. The absorber used had a mass ratio of MR=6.25. The MATLAB 
programs for modeling and analyzing this system are located in directory OFFCENTER in 
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Appendix A. Table 4-1 below compares absorber effectiveness for the compromise 
absorber location versus absorber locations optimized for a single mode: 





Absorber Located for Mode 3 
Effectiveness 


Absorber Located for Dual 
Mode Effectiveness 


Mode 3 Damping 
Ratio 


.00028 


.00023 


Mode 4 Damping 
Ratio 


0 


.00023 



Table 4-1: Dual Mode Absorber 0])timization 



The results indicate that relocation of the absorber causes some loss of mode 3 
damping ratio while raising mode 4 damping ratio from zero to a value equal to that of 
mode 3, This is the type of tradeoff would have to be performed in the design of a real 
marine cable absorber system. 

Figure 4-7 on the following page shows the mode 3 and 4 mode shapes with the 
absorber located at 7L/16. This figure demonstrates the utility of the STRING6 mode 
plotter in visualizing unorthodox mode shapes. 

So, absorbers can be located for effectiveness in both even and odd numbered 
modes. In choosing absorber location, the designer of a marine cable system would have 
to consider the excitation bandwidth and vibration modes within that bandwidth. If a 
single absorber location could not provide energy dissipation in all possible excitation 
modes, then multiple absorbers would need to be installed. 
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real(displacement mode shape) real(displacement mode shape) 



.q- 3 3rd MODE SHAPE. DAMPER AT 7L716 




-Iq- 3 4th MODE SHAPE, DAMPER AT 7L/16 




Figure 4-7: Mode 3 arul 4 Shapes with A!)sorl)er at 7L/16 
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4.5 Consideration of a Different Absorber Type 



Up to this point in the analysis, considerable insight has been gained in the 
behavior of the hinged, in-line, vibration absorber in the marine cable system. However, 
this type of absorber has not been able to produce the necessary damping ratio required to 
reduce vortex induced vibrations in a cable system in water. (See section 4.1). 

In this section, a simpler and more conventional absorber will be analyzed. This 
absorber is shown in figure 4-8 below: 




Figure 4-8: Han«iii«-Mass-Spriii‘' Daslipot Ahsoiber 
The absorber is of the conventional mass-spring-dashpot type. It is usually associated 
with discrete mechanical systems but could also be used to dissipate energy in a cable 
system. Li and Vandiver [10] derived the following transfer matrix for this absorber: 



[T] = 



1 

ry'm (j<yR + k) 
k - o)'m + jo)R 



0 

1 



(4-5) 
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A numerical model of a cable system using this absorber was constructed in 
MATLAB. The cable itself is identical to that of the baseline system of section 3.1. 
However, the absorber of figure 4-8 is located at the center of the cable in place of the in- 
line bending type absorber. The programs comprising this model are located in directory 
HANG in Appendix A. 

Den Hartog [1] formulated a method for choosing optimum absorber mass and 
stiffness for this type of absorber. His formulation was for discrete systems but modal 
analysis of the continuous cable system allows Den Hartog’ s method to be applied in an 
approximate way. The approach uses the following parameters: 
p = m / M = mass ratio = absorber mass/ main mass 
CO = Vk / m = natural frequency of absorber 

Q. = Vk / M = natural frequency of main system (cable without absorber) 

(0 / Cl - frequency ratio 

Den Hartog determined that the correct tuning of the mass spring dashpot 
absorber results from the frequency ratio being set as follows; 

f=l/(l+p) (4-6) 

This equation was used to set absorber stiffness. To do this, an example mass ratio of 
p = 0.2 was selected. The main mass was set at the modal mass of the cable which equals 
one half of total cable mass. D was set at the cable first mode natural frequency. 

Equation 4-6 then allowed computation of the absorber stiffness, k. 
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The following values were used in the analysis of the vibration absorber for the first 
vibration mode of the cable system: 
m = 5 kg = absorber mass 

M = 25 kg = cable modal mass = '/i of total cable mass 
|i = m / M = 0.2 

fi = 4.44 rad/sec = cable first mode natural frequency without absorber 
Using equation 4-6: 

f= 0.833 =m/0 



(0 = 3.70 rad/sec = 
k = 68.5 N/m, the correct stiffness for optimum tuning 

Once stiffness and mass were set, MULRUN was used to determine the optimum 
damping constant for the absorber designed in the procedure above. Figure 
4-9 below shows the resulting damping ratio performance and first mode shape: 




MODAL DAMPING RATIO VS. DAMPING CONSTANT 




Figure 4-9: Dam|)crO))timization for First Mode with Hanging Mass Spring Dashpot Absorber 
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It can be seen that optimum damping ratio is approximately 0.22, far higher than 
any damping ratio achieved with the in line absorber. This damping ratio is also above the 
critical value of = .04 necessary to reduce vortex induced vibration amplitude in water. 
This result also provides a reasonable reconciliation of results between the Den Hartog 
method and the transfer matrix method used in this thesis. Den Hartog states that the 
optimum damping constant for the mass spring dashpot absorber can be determined by the 
following formula; 

Rnptimum t 

= [3n/(8(l+p/)]'- (4-7) 

2 £2 m 

For a mass ratio of p = 0.2, this results in an R«piinu.mOf 9.25 N/(m/s). This corresponds 
roughly with the optimum damping constant of 15 N/(m/s) shown in figure 4-7. 

The hanging mass-spring-dashpot absorber represents a large improvement in 
absorber effectiveness over the bending type absorber discussed elsewhere in this thesis. 
This absorber has some possibility for use in marine cable systems and further analysis 
should be conducted. 
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Chapter 5 



Conclusion 

In this thesis, the transfer matrix method was used to construct various numerical 
models of marine cable systems. The performance of various vibration absorber 
configurations was studied. The following are the main findings of this research; 

1) Introduction of damping in the transfer matrix model causes the system 
characteristic equation to have complex roots. The use of a computer complex 
root solver was introduced and was shown to be useful in calculating system 
natural frequencies and damping ratios. 

2) Upon understanding the significance of complex natural frequencies, computer 
programs using transfer matrices were created to plot system mode shapes. 

The ability to visualize mode shapes was very important in understanding 
system behavior. 

3) The use of an in-line bending type vibration absorber was studied in some 
depth. A method was devised for tuning such an absorber to achieve 
maximum damping within a marine cable system. Variation of parameters such 
as absorber mass, length and number of absorbers was performed to study the 
effect of these variations on absorber performance. Even after optimization, 
the in-line absorber was unable to demonstrate sufficient damping to reduce 
vortex-induced vibrations in water. 
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4) The use of a hanging mass-spring-dashpot absorber was also studied for use in 
a marine cable system. This absorber was shown to have the ability to create 
sufficient damping to reduce vortex induced vibrations in water. 

The main focus of this thesis was to study the performance of the in-line bending 
type absorber in a finite length marine cable system. Although it was shown that this 
absorber provides insignificant damping to reduce vortex-induced vibrations in water, 
considerable insight was gained in this thesis. A working model was created which has the 
ability to calculate system natural frequencies, mode shapes and damping ratios. 
Additionally, it should be noted that the in-line absorber could be useful in reducing vortex 
induced vibrations in air. The studies in this thesis would be useful in optimizing the 
absorber for this application. 

Areas for potential research related to this thesis include: 

1) Use of the transfer matrix models in this thesis to study other marine cable 
system characteristics, such as the effects of distributed masses within the 
system. Li [5] performed considerable work in this area but the ability to 
visualize mode shapes would provide a good supplement his work. 

2) Extension of the analysis procedure in this thesis to other systems which can be 
represented by transfer matrices. 

3) Further study of the mass-spring-daspot absorber for use in marine cable 
systems. 
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Appendix A: MATLAB Files 



This appendix contains all MATLAB files used in this thesis. The files are divided 
into four directories, each corresponding to a different cable system model. The 
following is a list of these directories and a description of the system they model: 

CENTER; A cable system 50 meters long with an in line bending type absorber at 
its center. 

OFFCENTER; A cable system 50 meters long with an in line bending type 
absorber located at 7L/16. 

MULTABS; A cable system 50 meters long with three in line bending type 
absorbers, each located at a mode 3 antinode. 

HANG: A cable system 50 meters long with a hanging mass spring dashpot 
absorber located at its center. 

Variables and Dimensions 

The following variables are used throughout the MATLAB programs in this 
Appendix; 



Variable 


Corresponding Parameter 


Units 


a, b, c, d, e, f 


cable dimensions 


m 


Ml, M2 


absorber link masses 


kg 


T, H 


cable tension 


N 


K 


stiffness 


N/m 


R 


damping constant 


N/(m/s) 


11, 12 


absorber link moments of inertia 


kg m^ 


P 


cable mass density 


kg/m 


c 


cable wave propagation speed 


m/s 


k 


cable wave number 


m‘* 


w 


circular frequency 


rad/s 
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CENTER 

The following is a typical sequence for using the model CENTER to determine 
system natural frequencies and mode shapes and optimize system damping: 

1) In MATLAB, change directory to CENTER. 

2) Edit STRING 1 for desired absorber properties and STRTNG2 for desired cable 
properties. 

3) In STRING4 and STRINGl, ensure damping constant is set to zero. Run STRING4 
to plot ti 2 of [Toveraii] versus frequency. The zeroes of this plot are the undamped system 
natural frequencies. 

4) If desired, plot undamped system mode shapes as follows; 

a) For mode of interest, enter natural frequency in MATLAB window. 

b) Run STRING6 to plot mode shape. 

c) Repeat for other modes of interest. 

5) Perform absorber optimization for individual modes as follows; 

a) Edit MULRUN to iterate over desired range of damping constants. 

b) Also in MULRUN, set zrs to the undamped natural frequency for the mode of 
interest. This serves as a starting guess for the complex root solver. 

6) Set damping constant to global in STRINGl. Run MULRUN to plot real and 
imaginary natural frequency and damping ratio versus damping constant. The optimum 
damping constant is that which maximizes damping ratio. 

7) Plot damped system mode shape as follows; 

a) Set damping constant to the optimum value in STRINGl. 

b) Enter the complex natural frequency in the MATLAB window. 

c) Run STRING6 to plot mode shape. 

8) Repeat steps 5 through 7 for other modes of interest. 
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STRING 1 



% 

% This file contains the transfer matrix for a vibration 
% absorber. 

% Initial data 

T=5000; 

a=.03; 

b=.03; 

c=.2; 

d=.2; 

e=.2; 

f=.2; 

Ml=9.0; 

M2=M1; 

H=T; 

Il=.12; 

12 = 11 ; 

K=10; 

R=4000; 

h=c+d; 

l=e+f; 

al=(T*d)+(H*c)+(K*a^2)-(w'2*Il)+(i*w*R*b^2); 



a2=(T*f)+(H*e)+(K*a‘2)-(w'2*I2)+(i*w*R*b'2): 

a3=(K*a^2)+(i*w*R*b'2); 

gl=(-(Ml*l+M2*f)*al- 

(Ml*c*a3)+(w‘2*Ml*M2*f*c''2))/((M2*e*al)+(Ml*d+M2*h)*a3+(w‘2*Ml*M2*c*e*d)) 

g2=(-(M2*f*h^2)- 

(M 1 *l*d'2)+(a 1 *l+a3*h)/w'2)/((M2*e*a 1 )+(M 1 *d+M2*h)*a3+(w‘2*M 1 *M2*c*e*d)); 
g3=- 

((M 1 *l*d*c)+(al *l+a3*h)/w*2)/((M2*e*a 1 )+(M 1 *d+M2*h)*a3+(\v'2*M 1 *M2*c*e*d)); 

g4=-((M 1 *1+M2*f)*a3+(M 1 *c*a2)+(w^2*M 1 *M2*e*c*f))/((a2*h+a3*l)/w*2- 

(MPd*r2+M2*r2*h)); 

g5=((a2*h+a3*l)/w'2+(M2*e*f*h))/((a2*h+a3*l)/vv'2-(Ml*d*r2+M2*r2*h)); 

g6=(-a2*(M 1 *d+M2*h)-(M2*e*a3)+(w"2*M 1 *M2'-*=e*2*d))/((a2*h+a3*l)/w^2- 
(Ml*d*r2+M2*f2*h)); 
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t( 1 , l)=g 1 +g3=^((g4+g 1 *g6)/( l-g3*g6)); 
t( 1 ,2)=g2+g3 *((g5+g2*g6)/( 1 -g3 *g6)) ; 
t(2, 1 )=(g4+g 1 *g6)/( 1 -g3 *g6) ; 
t(2,2)=(g5+g2*g6)/(l-g3*g6); 



% 



STRING2 



% This file calculates the transfer matrix for a string of length L. 

%Initial data 

T=5000; 

L=25; 

p=1.0; 

c=sqrt(T/p); 

k=w/c; 

t(l,l)=cos(k*L); 

t(l,2)=sin(k*L)/(k*T); 

t(2,l)=-k*T*sin(k*L); 

t(2,2)=cos(k*L); 
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% 



STRING4 



% This file plots the overall T(l,2) for a system with 3 components. 
W=0; 

Ttotl2=0; 

R=0; 

for n= 1:250; 
w=.l*n; 

W(n)=w; 

string2 

Tl=t; 

T3=t; 
string 1 
T2=t; 

Ttot=Tl*T2*T3; 

Ttotl2(n)=Ttot(l,2); 

end 

plot(W,abs(Ttotl2)) 

grid 

xlabel(’w(rad/sec)’) 

ylabel(T(l,2)’) 



STRINGS 



% 



% This file contains the transfer matrix for a uniform bar. 
% Initial data: 

% bar length 
lb=.8; 

M=18.0; 

T=5000; 

Ib=M*lb*w'2/T; 

k=l/(2*(6+Ib)); 

t(U)=k*4*(3-Ib); 

t(l,2)=k*12*lbAT; 

t(2,l)=-k*M*w^2*(12-Ib); 

t(2,2)=k*4*(3-Ib); 
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STRING6 



% 



% THIS FILE PLOTS THE MODE SHAPES OF A STRING 

% WITH A DAMPER AT ITS CENTER. BEFORE RUNNING, 

% ENTER FREQUENCY IN THE MATLAB WINDOW AND 

% ENTER DAMPING IN STRING 1 . 

y=0; 

f=0; 

T=5000; 

P=1.0; 
c=sqrt(T/p); 
k=w/c; 
x=0:. 1:50.8; 



for n=l:25 1 

tl(l,l)=cos(k*x(n)); 

tl(l,2)=sin(k*x(n))/(k*T); 

tI(2,l)=-k*T*sin(k*x(n)); 

tl(2,2)=cos(k*x(n)); 

y(n)=tl(l,2); 

f(n)=tl(2,2); 

end 



string 1 

ts(l,l)=cos(25*k); 

ts(l,2)=sin(25*k)/(k*T); 

ts(2,2)=cos(25*k); 

ts(2,l)=-k*T*sin(25*k); 



for n=258;509 

tr(l,l)=cos(k*(x(n)-25.8)); 

tr(l,2)=sin(k*(x(n)-25.8))/(k*T); 

tr(2,l)=-k*T*sin(k*(x(n)-25.8)); 

tr(2,2)=cos(k*(x(n)-25.8)); 

tf=ts*t*tr; 

y(n)=tf(l,2); 

f(n)=tf(2,2); 

end 

plot(x(l;251),real(y(l;251)),x(258:509),real(y(258:509))) 

grid 

xlabel(’x location’) 
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- 2 - 



ylabel(’real(displacement mode shape)’) 
title(’3rd MODE MODE SHAPE’) 



5S 



% 



STRING 14 



% This file computes T(l,2) for a system with 3 components. 

function T12=stringl4(w) 

string2 

Tl=t; 

T3=t; 
string 1 
T2=t; 

Ttot=Tl*T2*T3; 

T12=Ttot(l,2); 
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% 



MULLER 



function [zros]=muller(fn,nprev,np,maxit,epl,ep2,fnreal,zros); 

% zeros=muller(fn,nprev,maxit,ep 1 ,ep2,fnreal,zeros) 

% This function is called by the driver program mulrun.m 
% Determines up to np zeros of the function specified by fn, 

% using quadratic interpolation, i.e., Muller’s Method. 

% This function is a translation of the Fortran routine given 
% by S. D. Conte and Carl de Boor, "Elementary Numerical 
% Analysis", McGraw Hill Book Company, 1980. The conversion 
% was made by Howard Wilson and Chris Roberts, Engineering 
% Mechanics Department, University of Alabama. 

% 

% To find more than one zero, function muller uses a procedure 
% known as deflation, which is implemented in function dflate. 

% fn: fuction that we want to find the roots. 

% zros: an array of initial estimates of all desired roots, set to zero if 
% no better estimates are available. 

% epl: relative error tolerance on the root. 

% epl: error tolerance on the function. 

% maxit: maximum number of iterations per root allowed. 

% np: number of roots desired. 

% nprev: number of roots previously computed, normally zero. 

% fnreal: a logical number; if TRUE, the program forces all approximations to 
% all the roots to be real. This makes it possible to use this routine 
% even if f(x) is defined only for real x. 

% Here we can take fnreal=l, if we only want to fine the real roots; we 

% can take fnreal=0 if we want to find both the real and complex roots. 

fn=’stringl4’; 

eps 1 =max(ep 1 , 1 -e- 1 2); eps2=max(ep2, 1 .e-20); 
for i=nprev+l:np; kount=0; 
while 1 

% Compute the first three estimates for zero as 
% zero, zero+0.5, zero-0.5 

zero=zros(i); h=0.00001+j *0.00001; zp=zero-fh; zm=zero-h; 

[fzr,dvdflp,kount,zros]=dflate(fn,zp,i,kount,zros); 

[fzr,fzrprv,kount,zros]=dflate(fn,zm,i,kount,zros); 

hprev=-2.0*h; dvdflp=(fzrprv-dvdflp)/hprev; 

[fzr,fzrdfl,kount,zros]=dflate(fn,zero,i,kount,zros); 

while 2 

% Do while kount<maxit or h is relatively big or fzr=f(zero) 
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% is not small or fzrdfl=fdeflated(zero) is not small or not 
% much bigger than its previous value fzrprv. 

divdfl=(fzrdfl-fzrprv)/h; divdf2=(divdfl-dvdflp)/(h+hprev); 
hprev=h; dvdflp=divdfl; c=divdfl+h*divdf2; 
sqr=c*c-4.0*fzrdfl*divdf2; 

if fnreal*real(sqr)<0.0, sqr=0.0; else, sqr=sqrt(sqr); end; 
if real(c)*real(sqr)+imag(c)*imag(sqr)<0, den=c-sqr; 
else, den=c+sqr; end 
if abs(den)<0.0, den=1.0; end 
h=-2.0*fzrdfl/den; fzrprv=fzrdfl; zero=zero+h; 
if kount>maxit, break, else 
while 3 

[fzr,fzrdfl,kount,zros]=dflate(fn,zero,i,kount,zros); 
if zros(i)==zero* 1.001, break; end 
% Check for convergence 

if abs(h)<epsl*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
% Check for divergence 

if abs(fzrdfl)>10.0*abs(fzrprv), h=h/2; zero=zero-h; 
else, break; end 
end; % while 3 
end; % end kount <= maxit 
if zros(i)==zero* 1.001, break; end 
if abs(h)<epsl*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
end; % while 2 
if kount>maxit, break; end 
if abs(h)<epsl*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
end; % while 1 
zros(i)=zero; 
end; % for loop 
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MULRUN 



% 



x=0; 
y=0; 
global R 

% input specs for Muller 
for n=l:51 
R=100*(n-1); 
x(n)=R; 

% tolerance, name of function that computes T12 etc: 



zrs = 121.; 
fn = ’string 14’; 
maxit = 500; 
epl = le-5 
% 

tolfcn = le-5; 
fnreal = 0; 
nprev = 0; 
np= 1; 



% Starting guess for the root, can be real or complex 
% name of program that computes T12 
% max # of iterations 

; %tolerance on the root (start 

values, le-2or so.) 

% relative tolerance on T12 
% = 0 complex search, = 1, real search 
% leave these as 0 
% and 1 



with 



ep2 = tolfcn*stringl4(zrs); % This sets the tolerance on T12 to be 

% le-5 of T12 estimated at your 
% guess 



% Call Muller 

root = muller(fn,nprev,np,maxit,epl,ep2,fnreal,zrs); 
y(n)=root; 
end 

zeta=imag(y)./real(y); 

subplot(31 1) 

plot(x,imag(y)) 

grid 

title(’4th mode’) 
xlabel(’R’) 

ylabel(’imag(nat freq)') 

subplot(312) 

plot(x,real(y)) 



lower 



starting 



62 



- 2 - 



grid 

xlabel(’R’) 

ylabel(’real(nat freq)’) 

subpIot(3l3) 

plot(x,zeta) 

grid 

xlabel(’R (N/(m/s))’) 
ylabel(’ modal damping ratio’) 
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DFLATE 



% 



% This is the subroutine needed by muller. 

function [fzero,fzrdfl,kount,zros]=dflate(F,zero,i,kount,zros); 

% (*fzero,fzrdfl,kount,zeros*)=dflate(F,zero,i,kount, zeros) 

% This function which is called by function muller is translated 
% from the Fortran routine given by S. D. Conte and Carl de Boor, 

% "Elementary Numerical Analysis", McGraw Hill Book Company, 1980. 
% The translation was made by Howard Wilson and Chris Roberts, 

% Engineering Mechanics Department, University of Alabama. 

kount=kount+l; fzero=feval(F,zero); fzrdfl=fzero; 

if i<2, return; end 

forj=2:i 

den=zero-zros(j-l); 
if abs(den)==0.0 
zros(i)=zero* 1.001; return 
else 

fzrdfl=fzrdfl/den; 

end 

end 
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OFFCENTER 



The following is a typical sequence for using the model OFFCENTER to 
determine system natural frequencies and mode shapes and optimize system damping; 

1) In MATLAB, change directory to OFFCENTER. 

2) Edit STRJNGIA for desired absorber properties and STRTNG2AL and STRING2AR 
for desired cable properties. 

3) In STRING4A and STRING 1 A, ensure damping constant is set to zero. Run 
STRING4A to plot t^of [Toveraii] versus frequency. The zeroes of this plot are the 
undamped system natural frequencies. 

4) If desired, plot undamped system mode shapes as follows; 

a) For mode of interest, enter natural frequency in MATLAB window. 

b) Run STRING6A to plot mode shape. 

c) Repeat for other modes of interest. 

5) Perform absorber optimization for individual modes as follows; 

a) Edit MULRUN to iterate over desired range of damping constants. 

b) Also in MULRUN, set zrs to the undamped natural frequency for the mode of 
interest. This serves as a starting guess for the complex root solver. 

6) Set damping constant to global in STRJNGIA. Run MULRUN to plot real and 
imaginary natural frequency and damping ratio versus damping constant. The optimum 
damping constant is that which maximizes damping ratio. 

7) Plot damped system mode shape as follows; 

a) Set damping constant to the optimum value in STRJNGIA. 

b) Enter the complex natural frequency in the MATLAB window. 

c) Run STRING6A to plot mode shape. 

8) Repeat steps 5 through 7 for other modes of interest. 



65 



I 



% 



STRING 1 A 



% This file contains the transfer matrix for a vibration 
% absorber. 



% Initial data 



T=5000; 

a=.03; 

b=.03; 

c=.2; 

d=.2; 

e=.2; 

f=.2; 

Ml=9; 

M2=iMl; 

H=5000; 

Il=.12; 

12 = 11 ; 

K=10; 

R= 100000; 

h=c+d; 

l=e+f; 



al=(T*d)+(H*c)+(K*a"2Hw"2*Il)+(i*w*R*b'2); 

a2=(T*f)+(H*e)+(K*a^2)-(w"2*I2)+(i*w*R*b^2); 

a3=(K*a'2)+(i*w*R*b'2); 



gl=(-(Ml*l+M2*f)*al- 

(M 1 *c*a3)+(w'2*M 1 *M2*f*c"2))/((M2*e*al)+(M 1 *d+M2*h)*a3+(w'2*M 1 *M2*c*e*d)) 
g2=(-(M2*f*h‘2)- 

(Ml*l*d"2)+(al*l+a3*h)/w^2)/((M2*e*al)+(Ml*d+M2*h)*a3+(w"2*Ml*M2*c*e*d)); 

g3=- 

((Ml*l*d*c)+(al*l+a3*h)/w^2)/((M2*e*al)+(Ml*d+M2*h)*a3+(w‘2*Ml*M2*c*e*d)); 
g4=-((M 1 *1+M2*f)*a3+(M 1 *c*a2)+(w'2*M 1 *M2*e*c*f))/((a2*h+a3*l)/vv"2- 

(Ml*d*r2+M2*f2*h)); 

g5=((a2*h+a3*l)/w"2+(M2*e*f*h))/((a2*h+a3*l)/w'2-(MPd*r2+M2*r2*h)); 
g6=(-a2*(M 1 *d+M2*h)-(M2*e*a3)+(w*2*M 1 *M2*e^2*d))/((a2*h+a3*l)/w'2- 

(Ml*d*r2+M2*f2*h)); 
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t( 1 , 1 )=g 1 +g3*((g4+g 1 *g6)/( 1 -g3 *g6)); 
t( 1 ,2)=g2+g3 *((g5+g2*g6)/( 1 -g3 *g6)) ; 
t(2, 1 )=(g4+g 1 *g6)/( 1 -g3*g6); 
t(2,2)=(g5+g2*g6)/(l-g3*g6); 
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STRING2AL 



% This file calculates the transfer matrix for a string of 
length L, 

% located on the left side of the absorber 

%Initial data 
T=5000; 

L= (7/16) *50; 

P=1 . 0; 

c=sqrt (T/p) ; 
k=w/c ; 

t(l,l)=cos(k*L) ; 
t(l,2)=sin(k*L) / (k*T) ; 
t (2 , 1) =-k*T*sin (k*L) ; 
t (2, 2) =cos (k*L) ; 
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% 



STRING2AR 



% This file calculates the transfer matrix for a string of length L, 
% located on the right side of the absorber. 

%Initial data 

T=5000; 

L=(9/16)*50; 

p=1.0; 

c=sqrt(T/p); 

k=w/c; 

t(U)=cos(k*L); 

t(l,2)=sin(k*L)/(k*T); 

t(2,l)=-k*T*sin(k*L); 

t(2,2)=cos(k*L); 
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% 



STRING4A 



% This file plots the overall T(l,2) for a system with 3 components. 
W=0; 

Ttotl2=0; 

R=0; 

for n= 1:200; 
w=.l*n; 

W(n)=w; 

string2al 

Tl=t; 

string2ar 

T3=t; 

string la 

T2=t; 

Ttot=Tl*T2*T3; 

Ttotl2(n)=Ttot(l,2); 

end 

plot(W,abs(Ttotl2)) 

grid 

xlabel(’w(rad/sec)’) 

ylabel(’T(l,2)’) 
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% 



STRING6A 



% THIS FILE PLOTS THE MODE SHAPES OF A STRING 
% WITH AN OFF CENTER DAMPER. ENTER FREQ. 

% FROM THE MATLAB WINDOW AND DAMPING 

% INSTRLNGIA 

x=0; 

y=0; 

T=5000; 

p=1.0; 

c=sqrt(T/p); 
k=w/c; 
x=0:. 1:50.8; 



for n= 1:220 

tl(l,l)=cos(k*x(n)); 

tl(l,2)=sin(k*x(n))/(k*T); 

tl(2,l)=-k*T*sin(k*x(n)); 

tl(2,2)=cos(k*x(n)); 

y(n)=tl(l,2); 

f(n)=tl(2,2); 

end 



string la 

ts(l,l)=cos(21.9*k); 
ts(l,2)=sin(21.9*k)/(k*T); 
ts(2,2)=cos(2 1 .9*k); 
ts(2,l)=-k*T*sin(21.9*k); 

for n=227:509 

tr(l,l)=cos(k*(x(n)-22.7)); 

tr( 1 ,2)=sin(k*(x(n)-22.7))/(k*T); 

tr(2,l)=-k*T*sin(k*(x(n)-22.7)); 

tr(2,2)=cos(k*(x(n)-22.7)); 

tf=ts*t*tr; 

y(n)=tf(L2); 

f(n)=tf(2,2); 

end 

plot(x( 1 ;220),real(y( 1 :220)),x(227 :509),real(y(227:509))) 
grid 

xlabel(’x location’) 
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ylabel(’real(displacement mode shape)’) 

title(’4TH MODE MODE SHAPE, DAMPER AT 3L/8’) 
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% 



STRING 14A 



% This file computes T(l,2) for a cable with an off center absorber. 

function T12=stringl4a(w) 

string2al 

Tl=t; 

string2ar 

T3=t; 

string la 

T2=t; 

Ttot=Tl*T2*T3; 

T12=Ttot(l,2); 






MULLER 



% 

function [zros]=muller(fn,nprev,np,maxit,ep 1 ,ep2,fnreal,zros); 

% zeros=muIIer(fn,nprev,maxit,epl ,ep2,fnreal, zeros) 

% This function is called by the driver program mulrun.m 
% Determines up to np zeros of the function specified by fn, 

% using quadratic interpolation, i.e., Muller’s Method. 

% This function is a translation of the Fortran routine given 
% by S. D. Conte and Carl de Boor, "Elementary Numerical 
% Analysis", McGraw Hill Book Company, 1980. The conversion 
% was made by Howard Wilson and Chris Roberts, Engineering 
% Mechanics Department, University of Alabama. 

% 

% To find more than one zero, function muller uses a procedure 
% known as deflation, which is implemented in function dflate. 

% fn: fuction that we want to find the roots. 

% zros: an array of initial estimates of all desired roots, set to zero if 
% no better estimates are available. 

% epl: relative error tolerance on the root. 

% ep2: error tolerance on the function. 

% maxit: maximum number of iterations per root allowed. 

% np: number of roots desired. 

% nprev: number of roots previously computed, normally zero. 

% fnreal: a logical number; if TRUE, the program forces all approximations to 
% all the roots to be real. This makes it possible to use this routine 
% even if f(x) is defined only for real x. 

% Here we can take fnreal=l, if we only want to fine the real roots; we 

% can take fnreal=0 if we want to find both the real and complex roots. 

fn=’ string 14a’; 

eps 1 =max(ep 1 , Le- 1 2) ; eps2=max(ep2, 1 .e-20) ; 
for i=nprev+l:np; kount=0; 
while 1 

% Compute the first three estimates for zero as 
% zero, zero+0.5, zero-0.5 

zero=zros(i); h=0.00001+j *0.00001; zp=zero+h; zm=zero-h; 

[fzr,dvdflp,kount,zros]=dflate(fn,zp,i,kount,zros); 

[fzr,fzrprv,kount,zros]=dflate(fn,zm,i,kount,zros); 

hprev=-2.0*h; dvdflp=(fzrprv-dvdflp)/hprev; 

[fzr,fzrdfl,kount,zros]=dflate(fn,zero,i,kount,zros); 

while 2 

% Do while kount<maxit or h is relatively big or fzr=f(zero) 
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% is not small or fzrdfl=fdeflated(zero) is not small or not 
% much bigger than its previous value fzrprv. 

divdfl=(fzrdfl-fzrprv)/h; divdf2=(divdfl-dvdflp)/(h+hprev); 
hprev=h; dvdflp=divdfl; c=divdfl+h*divdf2; 
sqr=c*c-4.0*fzrdfl*divdf2; 

if fnreal*real(sqr)<0.0, sqr=0.0; else, sqr=sqrt(sqr); end; 
if real(c)*real(sqr)+imag(c)*imag(sqr)<0, den=c-sqr; 
else, den=c+sqr; end 
if abs(den)<0.0, den=1.0; end 
h=-2.0*fzrdfl/den; fzrprv=fzrdfl; zero=zero+h; 
if kount>maxit, break, else 
while 3 

[fzr,fzrdfl,kount,zros]=dflate(fn,zero,i,kount,zros); 
if zros(i)==zero* 1.001, break; end 
% Check for convergence 

if abs(h)<epsl*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
% Check for divergence 

if abs(fzrdfl)>10.0*abs(fzrprv), h=h/2; zero=zero-h; 
else, break; end 
end; % while 3 
end; % end kount <= maxit 
if zros(i)==zero* 1.001, break; end 
if abs(h)<epsl*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
end; % while 2 
if kount>maxit, break; end 
if abs(h)<epsl*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
end; % while 1 
zros(i)=zero; 
end; % for loop 
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% 



MULRUN 



x=0; 
y=0; 
global R 

% input specs for Muller 
for n=l:21 
R=10000*(n-1); 
x(n)=R; 

% tolerance, name of function that computes T12 etc; 



zrs = 10.88; 
fn = ’string 14a’; 
maxit = 500; 
epl = le-5 

% 

tolfcn = le-5; 
fnreal = 0; 
nprev = 0; 
np = 1; 



% Starting guess for the root, can be real or complex 
% name of program that computes T12 
% max # of iterations 

; %tolerance on the root (start 

values, le-2or so.) 

% relative tolerance on T12 
% = 0 complex search, = 1, real search 
% leave these as 0 
% and 1 



with 



ep2 = tolfcn*stringl4a(zrs); % This sets the tolerance on T12 to be 

% le-5 of T12 estimated at your 
% guess 



% Call Muller 

root = muller(fn,nprev,np,maxit,epl,ep2,fnreal,zrs); 
y(n)=root; 
end 

zeta=imag(y)./real(y); 

subplot(31 1) 

plot(x,imag(y)) 

grid 

title(’4th mode’) 
xlabel(’R’) 

ylabel(’imag(nat freq)') 

subplot(312) 

plot(x,real(y)) 



lower 



starting 
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grid 

xlabel(’R’) 

ylabel(’real(nat freq)’) 

subplot(313) 

plot(x,zeta) 

grid 

xlabel(’R’) 

ylabel(’ modal damping ratio’) 



% 



DELATE 



% This is the subroutine needed by muller. 

function [fzero,fzrdfl,kount,zros]=dflate(F,zero,i,kount,zros); 

% (*fzero,fzrdfl,kount,zeros*)=dflate(F,zero,i,kount, zeros) 

% This function which is called by function muller is translated 
% from the Fortran routine given by S. D. Conte and Carl de Boor, 

% "Elementary Numerical Analysis", McGraw Hill Book Company, 1980. 
% The translation was made by Howard Wilson and Chris Roberts, 

% Engineering Mechanics Department, University of Alabama. 

kount=kount+l; fzero=feval(F,zero); fzrdfl=fzero; 
if i<2, return; end 
for j=2:i 

den=zero-zros(j- 1 ); 
if abs(den)==0.0 
zros(i)=zero* 1.001; return 
else 

fzrdfl=fzrdfl/den; 

end 

end 
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MULTABS 



The following is a typical sequence for using the model MULTABS to determine 
system natural frequencies and mode shapes and optimize system damping: 



1) In MATLAB, change directory to MULTABS. 

2) Edit STRJNGl A for desired absorber properties and STRING2AI and STRING2AO 
for desired cable properties. 

3) In STRJNG4A and STRINGIA, ensure damping constant is set to zero. Run 
STRING4A to plot tn of [Toveraii] versus frequency. The zeroes of this plot are the 
undamped system natural frequencies. 

4) If desired, plot undamped system mode shapes as follows; 

a) For mode of interest, enter natural frequency in MATLAB window. 

b) Run STRING6A to plot mode shape. 

c) Repeat for other modes of interest. 

5) Perform absorber optimization for individual modes as follows; 

a) Edit MULRUN to iterate over desired range of damping constants. 

b) Also in IVIULRUN, set zrs to the undamped natural frequency for the mode of 
interest. This serves as a starting guess for the complex root solver. 

6) Set damping constant to global in STRINGIA. Run N'lULRUN to plot real and 
imaginary natural frequency and damping ratio versus damping constant. The optimum 
damping constant is that which maximizes damping ratio. 

7) Plot damped system mode shape as follows; 

a) Set damping constant to the optimum value in STRINGIA. 

b) Enter the complex natural frequency in the MATLAB window. 

c) Run STRING6A to plot mode shape. 

8) Repeat steps 5 through 7 for other modes of interest. 
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% 



STRINGIA 



% This file contains the transfer matrix for a vibration 
% absorber. 

% Initial data 

T=5000; 

a=.03; 

b=.03; 

c=.2; 

d=.2; 

e=.2; 

f=.2; 

Ml=9; 

M2=M1; 

H=5000; 

Il=.12; 

12 = 11 ; 

K=10; 
global R; 
h=c+d; 
l=e+f; 

al=(T*d)+(H*c)+(K*a'2)-(w'2*Il)+(i*w*R*b^2); 

a2=(T*f)+(H*e)+(K*a^2)-(w'2*I2)+(i*w*R*b^2); 

a3=(K*a^2)+(i*w*R*b'2): 

gl=(-(Ml*l+M2*f)*al- 

(Ml*c*a3)+(w*2*Ml*M2*f*c*2))/((M2*e*al)+(Ml*d+M2*h)*a3+(w'2*Ml*M2*c*e*d)); 

g2=(-(M2*f*h'2)- 

(M1 *l*d*2)+(al *l+a3*h)/w‘2)/((M2*e*al)+(M 1 *d+M2*h)*a3+(w‘2*M 1 *M2*c*e*d)); 
g3=- 

((M 1 *l*d*c)+(a 1 *l+a3*h)/w‘2)/((M2*e*a 1 )+(M 1 *d+M2*h)*a3+(w*2*M 1 *M2*c*e*d)); 
g4=-((M 1 *1+M2*f)*a3+(M 1 *c*a2)+(w^2*M 1 *M2*e*c*f))/((a2*h+a3*l)/w*2- 

(Ml*d*r2+M2*r2*h)); 

g5=((a2*h+a3*l)/w‘2+(M2*e*f*h))/((a2*h+a3*l)/w‘2-(Ml*d*r2+M2*f2*h)); 
g6=(-a2*(M 1 *d+M2*h)-(M2*e*a3)+(w^2*M 1 *M2*e‘2*d))/((a2*h+a3*l)/vv*2- 

(Ml*d*r2+M2*r2*h)); 
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t( U )=g 1 +g3 *((g4+g 1 *g6)/( 1 -g3 *g6)) ; 
t( 1 ,2)=g2+g3*((g5+g2*g6)/( l-g3--^g6)); 
t(2, 1 )=(g4+g 1 *g6)/( 1 -g3 *g6); 
t(2,2)=(g5+g2*g6)/(l-g3*g6); 
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STRING2AI 



% 



% This file calculates the transfer matrix for a string of length L/3, 
% located on the left side of the abosorber 

%Initial data 
T=5000; 

L=(l/3)*50; 

p=1.0; 

c=sqrt(T/p); 

k=w/c; 

t(U)=cos(k=^L); 

t(l,2)=sin(k*L)/(k*T); 

t(2,l)=-k*T*sin(k*L); 

t(2,2)=cos(k*L); 
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% 



STRING2A0 



% This file calculates the transfer matrix for a string of length L/6, 
% located on the left side of the abosorber 

%Initial data 
T=5000; 

L=(l/6)*50; 

p=1.0; 

c=sqrt(T/p); 

k=w/c; 

t(l,l)=cos(k*L); 

t(l,2)=sin(k*L)/(k*T); 

t(2,l)=-k*T*sin(k*L); 

t(2,2)=cos(k*L); 
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% 



STRING4A 



% This file plots the overall T(l,2) for a system with absorbers located 
% at the mode 3 antinodes. Before running, set R to the desired value 
% in STRING 1 A. 

W=0; 

Ttotl2=0; 

R=0; 

for n= 1:100; 
w=100+.l*n; 

W(n)=w; 

string2ao 

Tl=t; 

string2ai 

T3=t; 

string la 

T2=t; 

Ttot=T 1 *T2*T3*T2*T3 *T2*T 1 ; 

Ttotl2(n)=Ttot(l,2); 

end 

plot(W,abs(Ttot 1 2)) 
grid 

xlabel(’w(rad/sec)’) 

ylabel(’T(l,2)’) 
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STRING6A 



% 



% THIS HLE PLOTS THE MODE SHAPES OF A STRING 
% WITH ABSORBERS AT THE MODE 3 ANTINODES. 

% ENTER FREQUENCY 

% FROM THE MATLAB WINDOW AND DAMPING 

% IN STRING 1 A 

x=0; 

y=0; 

T=5000; 

p=1.0; 

c=sqrt(T/p); 
k=w/c; 
x=0:. 1:52.4; 



for n=l:84 

tl(l,l)=cos(k*x(n)); 

tl(l,2)=sin(k*x(n))/(k*T); 

tl(2,l)=-k*T*sin(k*x(n)); 

tl(2,2)=cos(k*x(n)); 

y(n)=tl(l,2); 

f(n)=tl(2,2); 

end 

string la 



tso(l,l)=cos(8.33*k); 
tso( 1 ,2)=sin(8.33*k)/(k*T); 
tso(2,2)=cos(8.33*k); 
tso(2,l)=-k*T*sin(8.33*k); 

tsi(l,l)=cos(16.67*k); 

tsi(l,2)=sin(16.67*k)/(k*T); 

tsi(2,2)=cos(16.67*k); 

tsi(2,l)=-k*T*sin(16.67*k); 

for n=92:259 

tlc( 1 , 1 )=cos(k*(x(n)-9. 1 )); 

tlc( 1 ,2)=sin(k*(x(n)-9. l))/(k*T); 

tlc(2, 1 )=-k*T*sin(k*(x(n)-9. 1 )); 

tlc(2,2)=cos(k*(x(n)-9. 1)); 

tf=tso*t*tlc; 

y(n)=tf(l,2); 
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f(n)=tf(2,2); 

end 



for n=267;434 

trc( 1 , 1 )=cos(k*(x(n)-26.6)); 

trc( 1 ,2)=sin(k*(x(n)-26.6))/(k*T); 

trc(2, 1 )=-k*T*sin(k*(x(n)-26.6)); 

trc(2,2)=cos(k*(x(n)-26.6)); 

tf=tso*t*tsi*t*trc; 

y(n)=tf(l,2); 

f(n)=tf(2,2); 

end 

for n=442:525 

tr(U)=cos(k*(x(n)-44.1)); 

tr( 1 ,2)=sin(k*(x(n)-44. 1 ))/(k*T); 

tr(2,l)=-k*T*sin(k*(x(n)-44.1)); 

tr(2,2)=cos(k*(x(n)-44.1)): 

tf=tso*t*tsi*t*tsi*t*tr; 

y(n)=tf(l,2); 

f(n)=tf(2,2); 

end 

plot(x( 1 :84),real(y( 1 :84)),x(92:259),real(y(92:259)),x(267;434),real(y(267;434)),x(442:525),real(y(442:5 
grid 

xlabel(’x location’) 

ylabel(’real(displacement mode shape)’) 

title(’3rd MODE MODE SHAPE, DAMPERS AT ANTINODES’) 
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% 



STRING 14A 



% This file computes T(l,2) for a cable with absorbers located at the 
% mode 3 antinodes. 

function T12=stringl4a(w) 

string2ao 

Tl=t; 

string2ai 

T3=t; 

string la 

T2=t; 

Ttot=Tl *T2*T3*T2*T3*T^*T1 * 

T12=Ttot(l,2); 
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% 



MULLER 



function [zros]=muller(fn,nprev,np,maxit,epl ,ep2,fnreal,zros); 

% zeros=muller(fn,nprev,maxit,ep 1 ,ep2,fnreal, zeros) 

% This function is called by the driver program mulrun.m 
% Determines up to np zeros of the function specified by fn, 

% using quadratic interpolation, i.e., Muller’s Method. 

% This function is a translation of the Fortran routine given 
% by S. D. Conte and Carl de Boor, "Elementary Numerical 
% Analysis", McGraw Hill Book Company, 1980. The conversion 
% was made by Howard Wilson and Chris Roberts, Engineering 
% Mechanics Department, University of Alabama. 

% 

% To find more than one zero, function muller uses a procedure 
% known as deflation, which is implemented in function dflate. 

% fn: fuction that we want to find the roots. 

% zros: an array of initial estimates of all desired roots, set to zero if 
% no better estimates are available. 

% epl; relative error tolerance on the root. 

% ep2: error tolerance on the function. 

% maxit: maximum number of iterations per root allowed. 

% np: number of roots desired. 

% nprev: number of roots previously computed, normally zero. 

% fnreal: a logical number; if TRUE, the program forces all approximations to 
% all the roots to be real. This makes it possible to use this routine 
% even if f(x) is defined only for real x. 

% Here we can take fnreal=l, if we only want to fine the real roots; we 

% can take fnreal=0 if we want to find both the real and complex roots. 

fn=’ string 14a’; 

eps 1 =max(ep 1 , 1 .e- 1 2); eps2=max(ep2, 1 .e-20); 
for i=nprev+l:np; kount=0; 
while 1 

% Compute the first three estimates for zero as 
% zero, zero+0.5, zero-0.5 

zero=zros(i); h=0.00001+j *0.00001; zp=zero+h; zm=zero-h; 

[fzr,dvdflp,kount,zros]=dflate(fn,zp,i,kount,zros); 

[fzr,fzrprv,kount,zros]=dflate(fn,zm,i,kount,zros); 

hpre v=-2.0*h; d vdf 1 p=(fzrprv-d vdf 1 p)/hprev; 

[fzr,fzrdfl,kount,zros]=dflate(fn,zero,i,kount,zros); 

while 2 

% Do while kount<maxit or h is relatively big or fzr=f(zero) 
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% is not small or fzrdfl=fdeflated(zero) is not small or not 
% much bigger than its previous value fzrprv. 

divdfl=(fzrdfl-fzrprv)/h; divdf2=(divdfl-dvdflp)/(h+hprev); 
hprev=h; dvdflp=divdfl; c=divdfl+h*divdf2; 
sqr=c*c-4.0*fzrdfl*divdf2; 

if fnreal*real(sqr)<0.0, sqr=0.0; else, sqr=sqrt(sqr); end; 
if real(c)*real(sqr)+imag(c)*imag(sqr)<0, den=c-sqr; 
else, den=c+sqr; end 
if abs(den)<0.0, den=1.0; end 
h=-2.0*fzrdfl/den; fzrprv=fzrdfl; zero=zero+h; 
if kount>maxit, break, else 
while 3 

[fzr ,fzrdfl ,kount,zros]=dfl ate (fn ,zero,i ,kount,zros) ; 
if zros(i)==zero* 1.001, break; end 
% Check for convergence 

if abs(h)<epsl*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
% Check for divergence 

if abs(fzrdfl)>10.0*abs(fzrprv), h=h/2; zero=zero-h; 
else, break; end 
end; % while 3 
end; % end kount <= maxit 
if zros(i)==zero* 1.001, break; end 
if abs(h)<epsl*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
end; % while 2 
if kount>maxit, break; end 
if abs(h)<epsl*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
end; % while 1 
zros(i)=zero; 
end; % for loop 
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% 



MULRUN 



x=0; 
y=0; 
global R 

% input specs for Muller 
for n=l:21 
R=500*(n-1); 
x(n)=R; 

% tolerance, name of function that computes T12 etc: 



zrs = 104; 
fn = ’stringl4a’; 
maxit = 500; 
epl = le-5 
% 

tolfcn = le-5; 
fnreal = 0; 
nprev = 0; 
np= 1; 



% Starting guess for the root, can be real or complex 
% name of program that computes T12 
% max # of iterations 

; %tolerance on the root (start 

values, le-2or so.) 

% relative tolerance on T12 
% = 0 complex search, = 1 , real search 
% leave these as 0 
% and 1 



with 



lower 



ep2 = tolfcn*stringl4a(zrs); % This sets the tolerance on T12 to be 

% le-5 of T12 estimated at your starting 
% guess 



% Call Muller 

root = muller(fn,nprev,np,maxit,epl,ep2,fnreal,zrs); 
y(n)=root; 
end 

zeta=imag(y)./real(y); 

subplot(311) 

plot(x,imag(y)) 

grid 

title(’ mode, dampers at mode 3 antinodes’) 
xlabel(’R’) 

ylabel(’imag(nat freq)’) 

subplot(312) 

plot(x,real(y)) 
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grid 

xlabel(’R’) 

ylabel(’real(nat freq)’) 

subplot(313) 

plot(x,zeta) 

grid 

xlabel(’R’) 

ylabelC modal damping ratio’) 
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DFLATE 



% 

% This is the subroutine needed by muller. 

function [fzero,fzrdfl,kount,zros]=dflate(F,zero,i,kount,zros); 

% (*fzero,fzrdfl,kount,zeros*)=dflate(F,zero,i,kount, zeros) 

% This function which is called by function muller is translated 
% from the Fortran routine given by S. D. Conte and Carl de Boor, 

% "Elementary Numerical Analysis", McGraw Hill Book Company, 1980. 
% The translation was made by Howard Wilson and Chris Roberts, 

% Engineering Mechanics Department, University of Alabama. 

kount=kount+l; fzero=feval(F,zero); fzrdfl=fzero; 

if i<2, return; end 

forj=2:i 

den=zero-zros(j - 1 ); 
if abs(den)==0.0 
zros(i)=zero* 1.001; return 
else 

fzrdfl=fzrdfl/den; 

end 

end 
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HANG 



The following is a typical sequence for using the model HANG to determine 
system natural frequencies and mode shapes and optimize system damping: 

1) In MATLAB, change directory to HANG. 

2) Edit STRING! for desired absorber properties and STRING2 for desired cable 
properties. 

3) In STRING4 and STRING!, ensure damping constant is set to zero. Run STRING4 
to plot ti 2 of [Toveraii] versus frequency. The zeroes of this plot are the undamped system 
natural frequencies. 

4) If desired, plot undamped system mode shapes as follows; 

a) For mode of interest, enter natural frequency in MATLAB window. 

b) Run STRING6 to plot mode shape. 

c) Repeat for other modes of interest. 

5) Perform absorber optimization for individual modes as follows: 

a) Edit MULRUN to iterate over desired range of damping constants. 

b) Also in MULRUN, set zrs to the undamped natural frequency for the mode of 
interest. This serves as a starting guess for the complex root solver. 

6) Set damping constant to global in STRING!. Run MULRUN to plot real and 
imaginary natural frequency and damping ratio versus damping constant. The optimum 
damping constant is that which maximizes damping ratio. 

7) Plot damped system mode shape as follows; 

a) Set damping constant to the optimum value in STRING!. 

b) Enter the complex natural frequency in the MATLAB window. 

c) Run STRING6 to plot mode shape. 

8) Repeat steps 5 through 7 for other modes of interest. 
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% 



STRING 1 



% This file contains the transfer matrix for a hanging mass spring 
% dashpot vibration absorber. 

% Initial data 

M=5; 

K=68.45; 

R=15; 



t(l,2)=0; 

t(2,l)=-(w"2*M*(i*w*R+K))/(K-w"2*M+i*w*R); 

t(2,2)=l; 
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% 



STRING2 



% This file calculates the transfer matrix for a string of length L. 

%Initial data 

T=5000; 

L=25; 

P=1.0; 

c=sqrt(T/p); 

k=w/c; 

t(l,l)=cos(k*L); 
t(l,2)=sin(k*L)/(k*T); 
t(2, 1 )=-k*T*sin(k*L); 
t(2,2)=cos(k*L); 
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% STRING4 

% This file plots the overall T(l,2) for a system with 3 components. 
W=0; 

Ttotl2=0; 

R=0; 

for n= 1:200; 
w=0.1*n; 

W(n)=w; 

string2 

Tl=t; 

T3=t; 
string 1 
T2=t; 

Ttot=Tl*T2*T3; 

Ttotl2(n)=Ttot(l,2); 

end 

plot(W,abs(Ttotl2)) 

grid 

xlabel(’w(rad/sec)’) 

ylabel(’T(l,2)’) 
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% STRING6 

% THIS HLE PLOTS THE MODE SHAPES OF A STRING 
% WITH A MASS SPRING DASHPOT DAMPER AT ITS % 
BEFORE RUNNING, 

% ENTER FREQUENCY IN THE MATLAB WINDOW AND 

% ENTER DAMPING IN STRING 1 . 

y=0; 

f=0; 

T=5000; 

p=1.0; 

c=sqrt(T/p); 

k=w/c; 

x=0:.l;50; 



for n= 1:251 

tl(l,l)=cos(k*x(n)); 

tl(l,2)=sin(k*x(n))/(k*T); 

tl(2, 1 )=-k*T*sin(k*x(n)) ; 

tl(2,2)=cos(k*x(n)); 

y(n)=tl(l,2); 

f(n)=tl(2,2); 

end 



string 1 

ts(l,l)=cos(25*k); 

ts(l,2)=sin(25*k)/(k*T); 

ts(2,2)=cos(25*k); 

ts(2,l)=-k*T*sin(25*k); 



for n=252:501 

tr(l,l)=cos(k*(x(n)-25)): 

tr( 1 ,2)=sin(k*(x(n)-25))/(k*T); 

tr(2,l)=-k*T*sin(k*(x(n)-25)); 

tr(2,2)=cos(k*(x(n)-25)); 

tf=ts*t*tr; 

y(n)=tf(l,2); 

f(n)=tf(2,2); 

end 

plot(x(l:251),real(y(l:251)),x(252:501),real(y(252:501))) 

grid 



CENTER. 
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. T 



xlabel(’x location’) 

ylabel(’real(displacement mode shape)’) 
title(’3rd MODE MODE SHAPE’) 
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% 



STRING 14 



% This file computes T(l,2) for a system with 3 components. 

function T12=stringl4(w) 

string2 

Tl=t; 

T3=t; 
string 1 
T2=t; 

Ttot=Tl*T2*T3; 

T12=Ttot(l,2); 
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% 



MULLER 



function [zros]=muIIer(fn,nprev,np,maxit,epl,ep2,fnreal,zros); 

% zeros=muller(fn,nprev,maxit,epl,ep2,fnreal,zeros) 

% This function is called by the driver program mulrun.m 
% Determines up to np zeros of the function specified by fn, 

% using quadratic interpolation, i.e., Muller’s Method. 

% This function is a translation of the Fortran routine given 
% by S. D. Conte and Carl de Boor, "Elementary Numerical 
% Analysis", McGraw Hill Book Company, 1980. The conversion 
% was made by Howard Wilson and Chris Roberts, Engineering 
% Mechanics Department, University of Alabama. 

% 

% To find more than one zero, function muller uses a procedure 
% known as deflation, which is implemented in function dflate. 

% fn: fuction that we want to find the roots. 

% zros: an array of initial estimates of all desired roots, set to zero if 
% no better estimates are available. 

% epl: relative error tolerance on the root. 

% ep2: error tolerance on the function. 

% maxit; maximum number of iterations per root allowed. 

% np: number of roots desired. 

% nprev: number of roots previously computed, normally zero. 

% fnreal: a logical number; if TRUE, the program forces all approximations to 
% all the roots to be real. This makes it possible to use this routine 
% even if f(x) is defined only for real x. 

% Here we can take fnreal=l, if we only want to fine the real roots; we 

% can take fnreal=0 if we want to find both the real and complex roots. 

fn=’stringl4’; 

eps 1 =max(ep 1 , 1 .e- 12); eps2=max(ep2, 1 .e-20); 
for i=nprev+l:np; kount=0; 
while 1 

% Compute the first three estimates for zero as 
% zero, zero+0.5, zero-0.5 

zero=zros(i); h=0.00001+j *0.00001; zp=zero+h; zm=zero-h; 

[fzr,dvdflp,kount,zros]=dflate(fn,zp,i,kount,zros); 

[fzr,fzrprv,kount,zros]=dflate(fn,zm,i,kount,zros); 

hprev=-2.0*h; dvdflp=(fzrprv-dvdflp)/hprev; 

[fzr,fzrdfl,kount,zros]=dflate(fn,zero,i,kount,zros); 

while 2 

% Do while kount<maxit or h is relatively big or fzr=f(zero) 
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% is not small or fzrdfl=fdeflated(zero) is not small or not 
% much bigger than its previous value fzrprv. 

divdfl=(fzrdfl-fzrprv)/h; divdf2=(divdfl-dvdflp)/(h+hprev); 
hprev=h; dvdflp=divdfl; c=divdfl+h*divdf2; 
sqr=c*c-4.0*fzrdfl*divdf2; 

if fnreal*real(sqr)<0.0, sqr=0.0; else, sqr=sqrt(sqr); end; 
if real(c)*real(sqr)+imag(c)*imag(sqr)<0, den=c-sqr; 
else, den=c+sqr; end 
if abs(den)<0.0, den=1.0; end 
h=-2.0*fzrdfl/den; fzrprv=fzrdfl; zero=zero+h; 
if kount>maxit, break, else 
while 3 

[fzr,fzrdfl,kount,zros]=dflate(fn,zero,i,kount,zros); 
if zros(i)==zero* 1.001, break; end 
% Check for convergence 

if abs(h)<epsl*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
% Check for divergence 

if abs(fzrdfl)>10.0*abs(fzrprv), h=h/2; zero=zero-h; 
else, break; end 
end; % while 3 
end; % end kount <= maxit 
if zros(i)==zero* 1.001, break; end 
if abs(h)<epsl*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
end; % while 2 
if kount>maxit, break; end 
if abs(h)<epsl*abs(zero), break; end 
if max(abs(fzr),abs(fzrdfl))<eps2, break; end 
end; % while 1 
zros(i)=zero; 
end; % for loop 
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MULRUN 



% 



x=0; 
y=0; 
global R 

% input specs for Muller 

for n=l:41 

R=.5*n; 

x(n)=R; 

% tolerance, name of function that computes T12 etc: 



zrs = 4.07; % Starting guess for the root, can be real or complex 

fn = ’stringH’; % name of program that computes T12 
maxit = 500; % max # of iterations 

epl = le-5; %tolerance on the root (start with 

% values, le-2or so.) 
tolfcn = le-5; % relative tolerance on T12 

fnreal = 0; % = 0 complex search, = 1, real search 

nprev = 0; % leave these as 0 

np = 1 ; % and 1 



ep2 = tolfcn*stringl4(zrs); % This sets the tolerance on T 12 to be 

% le-5 of T12 estimated at your 
% guess 



% Call Muller 

root = muller(fn,nprev,np,maxit,epl,ep2,fnreal,zrs); 
y(n)=root; 
end 

zeta=imag(y)./real(y); 

subplot(31 1) 

plot(x,imag(y)) 

grid 

title(’4th mode’) 
xlabel('R’) 

ylabel(’imag(nat freq)’) 

subplot(312) 

plot(x,real(y)) 



lower 



starting 
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grid 

xlabel(’R’) 

ylabel(’real(nat freq)’) 

subplot(313) 

plot(x,zeta) 

grid 

xlabel(’R’) 

ylabelC modal damping ratio’) 
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% 



DFLATE 



% This is the subroutine needed by muller. 

function [fzero,fzrdfl,kount,zros]=dflate(F,zero,i,kount,zros); 

% (*fzero,fzrdfl,kount,zeros*)=dflate(F,zero,i,kount, zeros) 

% This function which is called by function muller is translated 
% from the Fortran routine given by S. D. Conte and Carl de Boor, 

% "Elementary Numerical Analysis", McGraw Hill Book Company, 1980. 
% The translation was made by Howard Wilson and Chris Roberts, 

% Engineering Mechanics Department, University of Alabama. 

kount=kount+l; fzero=feval(F,zero); fzrdfl=fzero; 
if i<2, return; end 
for j=2:i 

den=zero-zros(j- 1 ) ; 
if abs(den)==0.0 
zros(i)=zero* 1 .00 1 ; return 
else 

fzrdfl=fzrdfl/den; 

end 

end 
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Appendix B: Transfer Matrix Library 



The following appendix is a direct reproduction of the transfer matrix library 
derived by Li Li [5], It is included here so that this thesis can serve as a stand alone 
reference on transfer matrix analysis 
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Appendix A 

Transfer Matrix Library 



Definitions of the State Vectors and Transfer Matrices: (see Section 2.1) 

A.l Strings 

1) A string with varying tension 

The WKB solution to a string with varying-tension (see Fig. A-1) is given in [20] 
as 



y = 



1 



r[C:e^^ + C2e-^^] 






The transfer matrix can be found bom this solution as 



(A.l) 



Fi = 



ill ii2 
<21 <22 



(A.2) 



where: 






(p|r._,roi 



sinW 



ti2 = 






sinW 



'=■ = - Mp=Ti..roi + ji 



T'/ 'Tt 
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Figure A-1: A string with varying tension 






(<^T._,r.)T 



■sinW 



where W = w /g' — ^ — . If the tension varies linearly, then W 

V 

and T( = TU = T. 






2) A string with constant tension 

For a string with constant tension (see Fig. A-2), with k = \J ~ ^ the 
transfer matrix is 



Fi = 



coskle ^sinklc 
—kTsinklc coskl^ 



where 4 is the distance between the two state vectors at t — 1 and i. 



(A.3) 
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Figure A-2: A string with constant tension 




Figure A-3: A general, rigid lump 

A. 2 Rigid Lumps 

1) A general, rigid lump 

Section 2.5 gives the detailed derivation of the transfer matrix for a general, rigid 
lump (see Fig. A-3), the result is 




Ti{L, -l) + + (ic - O'] LI 

-Mu;^[Ti{L, -1) + Ti.^l - MwVO Ti{L, -l) + Ti.J - -f P] 

(A.4) 



where A = Ti{L, - /) + T._i/ - - {L, - 1)1]. 

2) Application 1: a uniform bar 

For a rigid, uniform bar {see Fig. A-4) with length li, and linear density pi, then 
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Figure A-4; A rigid, uniform bar 
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Figure A-5: A uniform, solid sphere 



Lc = h, I = M = pih, and 



The transfer matrix is 



Fi = 



ilii 

T 



4(3 - h) 

2(6 + 4) _ 4) 4(3 - 4) 



(A.5) 



where 4 = ^ average tension if there is a difference in tension at the 

two ends of the lump. 

3) Application 2: a solid sphere 

For a uniform, solid sphere (see Fig. A-5) with diameter D and volume density 
p„, then Le = D, I = ^, M — — and The transfer matrix is 



Fi = 



1 

20 -i 34 



20 - 74 

-Mo>2(20 - 24) 



200 

T 

20 - 74 



(A.6) 
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Figure A-6: A cone 

where I, = and T is the average tension if there is a difference in tension at 

two ends of the lump. 

4) Application 3: a cone 

For a cone (see Fig. A-6) with height h, bottom radius rj, and volume density 
then Lc = I = |/i, M = , and The transfer matrix is 



Fi = 



1 

1 - l,{rl - h^) 



1 - I,{rl Ih^) 

-Mu;^l - I^ri + i/i2) 



4h 

T.+3T._j 

1 - Ic(ri -I- 4h^) 



(A.7) 



where le = 



5(Ti+3Ti.i)h- 



5) Application 4: a point mass 

The transfer matrix for a point mass with mass M can be obtained directly from 
the transfer matrix for the general, rigid lump with = I = r = 0: 



Fi = 



1 0 
-Mu}^ 1 



(A.8) 
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Figure A-7: A general, rigid lump with damping devices (type I) 

A. 3 Damping Devices 

For a rigid lump with a damping device, the transfer matrix is defined the same 
way as that for a rigid lump without a damping device, 



/ \ 
y 

< 




in 


^12 


< 


f.-j 


i \ 


i 


^21 


<22 




Ty' J 



The expressions for the elements of the transfer matrix are as follows: 



(A.9) 



tn = 7 i +73 



74 + 7i7e 
1 - 7376 ’ 



tl2 = 72-T73 



75 + 7276 
1 - 7376 ’ 



^21 



74 + 7i76 
1 - 7378 ’ 



For the system shown in Fig. A-7, 7’s are given as follows: 



75 + 7276 
1 - 7376 

(A.IO) 



9 s { 9 i 9 i 2 ~ 93911) ~ 95(92911 + 9199) 
912(9296 — 9*99) + 99(9*99 + 53 ^ 6 ) 



(A.ll) 



72 



99(99 — 929 io ) ~ g8(g3gl0 + ^ 12 ) 
912(9296 — 9*99) ^9(54^5 + 9396 ) 



(A.12) 
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(A.13) 



73 



99(9397 — 9 s) ~i~ Pl2(g2g7 + 9 s) 
912(9296 ~ 949s) + 9s(949s + 5356) 



74 



51(55518 ~ 53513) + 517(5255 + 9 s 9 s) 

5is(5357 ~ 55) + 513(5257 + 9 s) + 514(5255 + 9 s 9 s) 



(A.14) 



75 



5s5i3 — 515(5255 + 9 s 9 s) — 5»5i8 

518(5357 ~ 9 s) + 513(5257 + 53) + 514(5255 + 5353) 



(A.15) 



513(5256 ~ 545s) — 5ie(5255 + 535s) + 513(5455 + 9 s 9 s) 
513(5357 - 55) + 513(5257 + 53) + 514(5255 + 5353) 

where ^’s are given by 



(A.16) 



g, = - ^), ^2 = + Mj^), ^3 = + Mz^), 

54 = u:nh(l - k), ^ 4. ^6 = ^ H- <^^Mzh(‘i - 1), 

57 = 58 = %) 59 = ^ + ^ 

5io = C) 5ii = ^ ~ w^iWi^2(l — (^)i 5i2 = 

5i3 = ^ ~ ^14 = /4, ^15 = /a, 

5i6 = <j^^MzU(l ~ "f ) + ^» 5i7 = ^ M\lz(\ — (^), 5i8 = ^ — ^ 

ai = T- Hilz + A’ldj — u)^Ii + j(j}R\b\ 

0-2 = Hrl^ + Hllz + AjCj + K2CL2 J<*’(-^1^1 -^2^2) ~ ^^-^2 

0:3 = T{Iq + Hrh + K2CL2 ~ ^^-^3 + Jt»li?2^2 
/?1 = Ajdj + jlj}Rib\, /?2 = K2<i\ + JUli?2^2 



where if/ and ff, are the forces along the structured length at the left and right 
joint, respectively, which can be found by statics. 

For the system shown in Fig. A-8, 7’s are given as follows: 



-(Ml/ + A/2/)ai - M^caz + 

Maeai + (Afid + M2h)az + w^MiMaced ’ 



_ -M2fk^ - Mil<P + (aj + Qzh)l^\ 
■}■ (M\d, -}- M2h)o.z A AI\ M2^^ 

(A.17) 
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7s, n = 9i9iz + 9z9[z ~ 9\i{929i + 9z9i) ~ 9i9\i ~ 9i9'iz) 

7 e,n = ^13(^256 ~ 9^99) + 513(^2^6 - gig's) — g'leigzgs + gzgs) - ^16(^255 + gzg's) 
+g'isigig 3 + gzge) + 518(5455 + 535 s) 

7 i,d = 5(3(5357 - 55) - gisg's + 5(3(5257 + 5s) + 5i35s + 514(525^ + 535s) 

7s,<i ~ 7e,(i ~ 7i,d 

5 ' = jubW + i), g'e = g's = j<^b‘^2^ 

g'9 = + 3 )» 5(1 = 512 = 

5(3 = + t)> 5(s = 5(7 = 

g'ls = + i) 

For the system shown in Fig. A- 8 , 7 ’s are given by 



, jW{[Mi(:+c)+Af5/hi,.i+[Afi<i+.V/5(e+K)hi.„} _ ;‘-4’{-^7j..i+[A/ii+iV5(<+K)]7j.n} 

7i = 7T. ’ 7: - ::2 



'>'1.4 



■J.4 






73 = 






i,d 



74 = 






7s j 7s 4^3 

~5,«i 



A. 4 Supports 



For the support shown in Fig. A-9, the transfer matrix is 



Fi = 



0 

z 

(L^-z)z Lc-i , 



5 r ^ 0, Xc 



(A.26) 



where Ig is the mass moment of inertia of the rigid lump about the pinned point. 
For the support shown in Fig. A-10, the transfer matrix is 



f.= 



Lr-Z 

Z 



IqU:^-T{Lc-z) 

(Lc-xjx Lc-z , 



, 2 ^ 0 , X<- 



(A.27) 



where is the mass moment of inertia of the rigid lump about the pinned point. 



113 




Figure A-9: An intermediate support 




Figure A-10: An end support 
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Appendix C: 



Analysis of a Single Degree of Freedom System Using the 
Transfer Matrix and Complex Natural Frequency Approach 



The System: 










o o 

The Transfer Matrix Relates the State Vectors at Points 1 and 2: 



X: 




'Til 


Ti:' 


XI 


F: 




J-' 


1 

r \ 
r X 

H 


_Fi_ 



Derive the Transfer Matrix of the System Between Points 1 and 2: 



-Fi = K(X2 - Xi) + joR(X2 - Xi) = (K + jcoR)(X2 - Xj) 
solving for X 2 : 

X 2 = -F|/(K + jcoR) + Xi (2) 

F, + F2=M(dVdr) 

Fi + F 2 = -C0‘MX2 
writing F 2 in terms of Fi and X\: 

¥2 = -F, - coNI(-Fi/(K + jcoR) + X,) 

F 2 = -(o'Mxi + (o)'M/(K +jcoR) -1)F, 
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( 3 ) 



Equations 2 and 3 result in the transfer matrix: 



Til Ti 
T:i T; 



1 -1/(k + jV4>R) 

— ty'M -l + ry'M/(K + jry r) 



Applying Boundary Conditions: 

Xi = 0 due to system being attached to wall 

p 2 = 0 because we want to study the free response of the system 

so, equation I becomes: 



1 

X 

1 -J 

1 




1 

H 


1 

r 1 


"0 ■ 


1 

o 

1 




_T:i 


T:: 


1 

1 



yields the characteristic equation: 

T22 = 0 

- 1 + ry ■ M / (k + j r) = 0 

-co‘M+jcoR + K = 0 (4) 

The Complex Frequency Root: 

applying quadratic formula \yith 4KM < R‘: 

V4KM-R' . R 

From Section 3.3.1 of thesis, a damped harmonic response can be written: 
x(t) = A cos(cOrcni t) e^'"’"’-^ ‘ 



x(t) = A cos 
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Comparison With Traditional Method: 



In the traditional method, the characteristic equation is the same as 
equation 4; 

-o)^M +j(oR + K = 0 

but the complex root is not solved for. Instead the following are defined: 



coo = 




R 

ImM 



the undamped natural frequency (7) 



, the damping ratio (8) 



(Hi = \ , the damped natural frequency (9) 

The displacement solution is then written as: 
x(t) = A cos( codt ) exp(- 4" ®ot) (10) 

When equations 7, 8 and 9 are used to write equation 10 in terms of 
K, R and M. The result is: 



x(t) = A cos 



V4KM-R- 



2M 



exp 



R 

'2M 



t (11) 



which is the same as equation 6, the result obtained using the complex 
frequency approach. 



) 
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